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Abstract 

We investigate the dynamics of many interacting bubbles in boiling water by using a laser 
scattering experiment. Specifically, we analyze the temporal variations of a laser intensity 
signal which passed through a sample of boiling water. Our empirical results indicate that the 
return interval distribution of the laser signal does not follow an exponential distribution; 
contrariwise, a heavy-tailed distribution has been found. Additionally, we compare the 
experimental results with those obtained from a minimalist phenomenological model, finding 
a good agreement. 

1. Introduction 

Bubbles are ubiquitous in nature and their dynamics is both fascinating and very complex 
It is not surprising that bubbles are an effervescent source of research. For instance, bub- 
bles appear in the context of energy generations, 4|, collapsing bubbles^ can emit ligthjsl, 
turbulent thermal convection has been observed in a single soap bubblejGi], cooperative and 
avalanches-like dynamics are present in collapsing of aqueous foams 7|, |8|, and singularities 
emerge when air bubbles detache from a nozzle submerged in water 9j. 

Despite the fact that the fundamental equations ruling the behavior of moving fluids 
are well known, an analytical or even a numerical approach can become infeasible for many 
common situations. This is particularly true for bubbles in turbulent fluids, where, at higher 
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Reynolds numbers, the number of mesh points required to solve each bubble as well as the 
flow around it grows up leading to a long simulation time 10 1. 

A very familiar case, where bubbles appears in such contexts, is the boiling process of 



water 



For exemple, when the temperature reaches 100°C the vapor pressure is 1 bar 



and we can observe the spontaneous process of bubble formation (nucleation) 12j . Although 
it is an ordinary process, nucleate boiling has several complex aspects involving thermal 
interactions between bubbles and the heated surface and among the nucleation sites. There 
are also hydrodynamics interactions bubble to bubble and bubble to liquid bulk[l3]. In 
this scenario, simple models, whether phenomenological or not, and simple experiments 
have been designed to try to clarify this intricate dynamics. For instance, extensive studies 
have been done by considering nonlinear models and experiments of boiling 13|] as well as 
evaporation in microchannels|l41] or in short capillary tubejlS]. However, as far the authors 
know, much less attention has been paid to clouds of bubbles, i.e., a system containing many 
interacting bubbles (see for instance Ref. 16j). Our main goal here is attempt to fill this 
hiatus by using a simple experiment that basically consists of a laser beam passing through 
a sample of boiling water. We also comfront the experimental results with a minimalist 
model towards improving our understanding of this complex system. 

This article is organized as follows. Primarily, we describe the experimental setup and 
the data acquisition. Next, we report a statistical analysis of the data and present a model. 
Finally, we end this work with some concluding comments and a summary. 



2. Experimental setup and data presentation 

The experiment consists of samples of approximately 300 ml of distilled water at atmo- 
spheric pressure and confined by impermeable metallic walls with glass windows through 
which the laser beam is transmitted, as shown in Figured] (left panel). At the bottom, the 
confining vessel is in contact with a heat plate that employs a power of around 300 W, in 
such way that the temperature in the contact interface is approximately 300 °C. After the 
boiling process becomes stable, i.e., the water temperature stabilizes, we start to record the 
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Figure 1: (left panel) A schematic representation of the experimental setup: a laser beam passes through 
a sample of boiling water which is in contact with the heat plate. The laser intensity is monitored by 
the detector (photodiode) and recorded, (right panel) A typical laser intensity signal for a sample of 
approximately 300 ml of boiling water. 

He-Ne laser (10 mW) intensity signal that passes through the sample by using a photodi- 
ode detector (Thorlabs DETIOOA) coupled to an oscilloscope (Tektronix TDS5032B) with 
a sampling rate of a thousand points per second. Typical recording times are of the order of 
10 minutes and the height of the incident beam does not significantly change the statistical 
results, avoiding the water-air interface. 

Figure [T] (right panel) shows a typical record signal. We can see that the signal is char- 
acterized by intermittent valleys. While the emerging dynamic is complex, the individual 
processes generating the bubble are qualitatively simple. At the button of the confining 
vessel the temperature is higher what makes a small fraction of the liquid to evaporate, pro- 
ducing the bubbles. These nucleation sites are not static and depend on the heat transfer 
and also on the liquid-wall interactions. The bubbles depart from the nucleation sites and 
rise through the confining vessel. Along this movement the bubbles continuously interact 
with each other, with the liquid and with the walls. When one or more bubbles crossing the 
laser path they scatter the light, producing a decrease in the record signal. Naturally, this 
intrincated signal reflects the complex collective dynamical behavior of the bubbles. Similar 
situations are customary when dealing with time series. For instance, the earth seismic and 
geomagnetic activity can be investigated by considering a seismogram and the DST index. 
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3. Data analysis 



Due the conditions of the signal, a natural variable to investigate the system dynamics 



is the time difference between the extreme events characterized by sharped val 
analysis is frequently employed in the physics and financial literature 
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20| and it 
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shows to be useful when investigating the underling mechanism ruling the system 
24|. 

A possible manner to obtain these extreme events is considering a threshold value q 
and archiving the initial times ti for which laser signal is below this edge. The difference 
between two consecutive times, Tj = tj+i — tj, is the so called return interval. This procedure 
is presented in Figure [T], where the right panel shows horizontal line segments representing 
the return intervals for q = 0.5. Figure [2^ displays the probability density functions (pdf) 
of Ti, p{t), for three values of q. Clearly, the distribution is dependent on q, and it is 
also well known that for Gaussian random uncorrelated variables the return interval Ti is 
exponentially distributed according to 
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(1) 



where fq is the mean value of return interval related to the threshold value q. Figure [2^ 
shows our data compared with this pdf, resulting in a poor agreement. 

The weak agreement indicate that memory effects can be present in the bubble dynamics. 
To address this question, we may use the detrended fluctuation analysis (DFA) 26|. It 
basically consists of calculating the root mean square fluctuation function F{n) (see for 
instance j^^]) for the integrated and detrended time series for different values of the time scale 



n. When we have scale- invariant time series, F{n) follows a power law behavior F{n) 



~ n 



where h measures the degree of correlation in the time series: if h = 0.5, the series is 
uncorrelated, while h > 0.5 indicates long-range correlations. Figure [2]d presents the results 
concerning the laser signal where we found h ^ 0.65 leading to long-range correlations. 
Empirical results have claimed that in the presence of power law correlation in the data 



the pdf p(r) is usually adjusted by a stretched exponential 
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23| or by a WeilbuU 




Figure 2: (a) Probability distribution of the return interval r for three values of q (indicated in the figure) 
in comparison with the exponential distribution of EqlT] for fo.i = 0.0808 (continues line), fo.e = 0.0564 
(dashed line) and fo.s = 0.0438 (dashed-dotted line), (b) DFA analysis for the laser intensity signal: 
log]^Q[_F(n)] versus \ogn^{n) in comparison with a linear fit, where we found F{ri) oc with h w 0.65. 
(c) Distribution of scaled variable ^ = r/fg for eight equally spaced threshold values from 0.1 to 0.9 in 
comparison with a stretched exponential (continuos line) and with a WeilbuU distribution (dashed line) 
both with 7 = 2(1 — /i) = 0.7. (d) Distribution of the return interval increments = — compared 
with EqOwhen considering the stretched exponential (continuos line) and the WeilbuU distribution (dashed 
line) both with 7 — 0.7. Notice that in the tail of the distributions the noise increases. This occurs because 
of the finite size of the set {ri} and also due to the small probability of finding large time intervals. 



distribution 



m, i.e., 

p(r) ~ e-^("/^')' or p(r) ~ (r/f^)^-' e-^^"/^")' , (2) 



where A and B are constants and 7 is the exponent of the power law autocorrelation func- 
tion. These two distributions also emerge in the analytical approach of Santhanam and 



n 



Kantz[25'] when considering a Gaussian fractional noise with autocorrelation exponent 7. It 
is interesting to note that by employing the scaled variable ^ = r/fq both distributions be- 
come independent on q and that the exponents 7 and h are related via 7 = 2(1 — h). Figure 
Wp- displays the distributions of the scaled variable ^ and also the stretched exponential and 
the Weilbull distributions, which, due to the normalization and the unitary mean value of 
^, have just one parameter, 7, determined from h. From this figure, we observe a good data 
collapse but a poor agreement with both previous distributions, specially for small ^. In 
this case, the distributions of Eql2] underestimate the data for ^ G [0, 2] and overestimate 
for ^ G [3, 7]. The agreement is not improved if we find 7 via least square method. 

Another aspect to be investigated is the scaled return interval increments = ^j+i — ^j. 
We know that the pdf of the difference between two independent random numbers X — Y 
is given by the cross-correlation 28 1 



where X is distributed according to fx{x) and Y according to /y(a;). In the case of A,^, 
both distributions are the same and we can use a stretched exponential or the Weibull of 
EqJ2]to compare with the data. Figure shows this comparison, again finding an imperfect 
agreement. In addition, we have to mention that the return interval series is week correlated 
{h ^ 0.55), thus the previous equation should be viewed as an approximation. 

4. Modeling 

In principle, we should be able to describe the dynamics of boiling fluid since the physical 
transport phenomena are well known to follow the Navier-Stokes equation. However, tech- 
nical difficulties as numerical instability because of the complex boundary of the two phases 
and the large amount of simulation time required to do the integrations make this task 
very difficult. Thus, our goal is to understand this complex phenomenon from a minimalist 
model. Therefore, we retain only the relevant ingredients to reproduce the main aspects of 
the experimental data. 




(3) 
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In an approximate scenery, we can employ a two-state model for which only matters 
whether there are bubbles in the laser path or not. When there are bubbles, the laser bean 
is considered totally scattered and the passing intensity signal is zero. On the other hand, 
when there are no bubbles, the beam passes without scattering and the intensity signal is 
one. 

In addition to the above two-state approximation for the laser signal, a significant in- 
gredient of these empirical data is the power law correlation presented in Figure It has 
been reported that a symbol series (here or 1) can present lon g-ra nge correlations when 



it is generated by using only two uncorrelated random numbers [29|, |30||. In this context, 
a possible way to model our data is considering a one- dimensional lattice where the sites 
represent the laser intensity signal in a given time. In the first site of the lattice we start 
to fill it drawing a discrete random number from a Bernoulli distribution with probability 
parameter p. If the random number is 1, there are bubbles in the laser patch and the laser 
signal is zero for this time. Otherwise, there are no bubbles in the laser path and the laser 
signal is one for this time and also for the next [x] times (sites), where [x] is the integer part 
of a random number x distributed according to a distribution P{x). All of the other sites of 
this lattice are filled by repeating the above procedure. Due to the fact that the simulated 
laser signal is only zeros or ones, it does not depend on the threshold value q. Actually, in 
this approach the return interval r is exactly the length of the consecutive zero sites. Thus, 
we effectively focus our attention on the the scaled variable ^. 

The proposed model appears to be very ad hoc since we have a function to choose plus a 
parameter. However, when performing the simulations, we empirically found that short tail 
distributions such as exponential, Gaussian, lognormal, and gamma are not able to improve 
the agreement found when comparing with the analytical expressions of EqlH In contrast, 
when considering the most simple distribution with power law tail, i.e., the Pareto one 

{ak°'x~°'~^ if a; > fc, 
(4) 
if a; < /c, 

where A; > 0, a > are parameters of the distribution, the agreement with our data is very 
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Figure 3: A comparison between the experimental (a) and the simulated (b) return intervals. The simulated 
result was obtained by considering a Pareto distribution with k ^ I and a = 1.8 and the parameter p = 0.2 
for a lattice with 5 x 10^ sites averaged over 200 realizations. Figure (c) shows the distribution of the 
return interval for eight equally spaced threshold values from 0.1 to 0.9 compared with the simulated result 
(continuos lines). Figure (d) presents the comparison for the return interval increments = — ^j. 



good. 

Figure [3] shows the simulated return interval in comparison with experimental data. For 
the simulation, we have fixed the value k = 1 and the size of lattice in a such way to have 
5 X 10^ return intervals, typically the number found in the experimental data. To obtain the 
best fit parameters, we incrementally update the values of a and p getting the distribution 
p(^) which is averaged over 200 realizations and so confronted with the experimental data 
via the method of least squares. The best values found for these parameters are a = 1.8 
and p = 0.2. 



8 



The previous model claims for an analytical approach. Indeed, it may be view as a sum 
of random numbers in which the number terms is also a random number. Thus, supposing 
known the distribution of Pareto sums, we may write 



p(r) = iV^(l-p)"5„K,A;, 



(5) 



n=l 



where is the normalization factor, (1 — p)" is the probability of summing n consecutive 
numbers and Sn{anikn) is the distribution of Pareto sums. However, to obtain a general 



expression for S'„(a„,/c„) may not be a easy task 



31 



33|. It can be argued that these 



cumbersome calculations are avoided if we consider a stable distribution for P(x). In fact, 
for this case the sum of stable variables is well known to be a stable distribution, but the 
support of P{x) is only positive when the stable index, a', is less that one and the skewness 



34| . The stable distribution is asymptotically a power law with 



parameter is equals to one 

P{x) ~ a;~^~"' which is not compatible the exponent found when considering the Pareto 
distribution. 



The comparison with the previous model also indicates that the laser signal is very close 
to a point possess, since the deflection times are very short when compared with inter-event 
times Tj. By comparing our model with Refs. 35, Q, Is^], we find that the laser signal can 
be viewed as non-Poisson renewal process. This process is characterized by a sequence of 
events spaced by time intervals that are independent random variables. Moreover, the time 
intervals are draw from the same probability density iP{t) ~ r^'^, where yU is the renew index. 
Confronting this expression with the Pareto distribution, we can see that the renew index 
is /i = a + 1 ~ 2.8 for our case. In addition, following Grigolini et. al 35|], it is also possible 
to show that a fully asymmetric Levy stable distribution of index 5 = — 1)"^ ~ 0.56 
emerges for number of events in an time interval r'. Further, the renew process may be also 
connected with a noise whose power spectrum is 1//'', where r/ = 3 — yU ~ 0.2 in our case. 

The previous findings suggest that a 1//°'^ noise is present in the context of a first-order 
phase transition, in contrast with the common association among inverse power laws noise 
and critical phenomena 38 



This has been also re por ted in the context of first-order 

electronic phase transitions j40| and for polymer folding 4l|. From a general point of view, 
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all the empirical discoveries suggest that, at boiling temperature, the fluctuations between 
the two phases (here, the bubble and the liquid bulk) play an essential role in the system 
dynamics, generating the nontrivial aspects reported here. 

5. Summary 

We have reported statistical analysis from the bubbles dynamics in a sample of boiling 
water. Our analysis was focused on an experiment of laser scattering in which a laser 
beam passed through the boiling fluid having its intensity monitored. By using this time 
series we evaluated the return interval distribution flnding a non-exponential distribution. 
In addition, we verifled that the dynamical processes that generates the bubbles introduces 
nontrivial correlations. Employing a minimalist phenomenological model we were able to 
reproduce the experimental behavior successfully . The model also seems to suggest that a 
fundamental ingredient generating the nontrivial dynamics is the power law tail related to 
the waiting time for the bubbles passing through the laser path and also the correlations 
introduced by the draw of the two random numbers. 
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